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Abstract 

This paper considers Schrodinger operators, and presents a prob- 
abilistic interpretation of the variation (or shape derivative) of the 
Dirichlet groundstate energy when the associated domain is perturbed. 
This interpretation relies on the distribution on the boundary of a 
stopped random process with Feynman-Kac weights. Practical com- 
putations require in addition the explicit approximation of the normal 
derivative of the groundstate on the boundary. We then propose to use 
this formulation in the case of the so-called fixed node approximation 
of fermion groundstates, defined by the bottom eigenelements of the 
Schrodinger operator of a fermionic system with Dirichlet conditions 
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on the nodes (the set of zeros) of an initially guessed skew-symmetric 
function. We show that shape derivatives of the fixed node energy 
vanishes if and only if either (i) the distribution on the nodes of the 
stopped random process is symmetric; or (ii) the nodes are exactly 
the zeros of a skew-symmetric eigenfunction of the operator. We pro- 
pose an approximation of the shape derivative of the fixed node energy 
that can be computed with a Monte-Carlo algorithm, which can be 
referred to as Nodal Monte-Carlo (NMC). The latter approximation of 
the shape derivative also vanishes if and only if either (i) or (ii) holds. 



1 Introduction and results 



Throughout this paper, we consider a Schrodinger operator in 
form: 

H = ~ + V, 



of the 
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with a smooth potential V going to infinity at infinity, and acting on real 
valued functions generically denoted with the letter if) ('wave functions'). 
Such functions ip will be defined up to a real valued multiplicative constant 
(e.g. in eigenvalue and/or variational problems). We also consider a general 
family 

e h-> n e (2) 

of open smooth domains in R d depending sufficiently smoothly of a parameter 
€ MP. The boundary will be denoted dVLg. Gradients in the space M. d will 
be denoted V, and gradients with respect to 9, Vg. The Dirichlet groundstate 
and its Dirichlet groundstate energy (ipg, Eq), are then defined as the unique 
bottom eigenelements of H, solution to the variational problem: 



En := inf 



V 



ifiH (iP) 



\ 



^\dn g = 



L 

Ja, 



*\2 



(3) 



Calculus of variations detailed in Section 2 then yields the shape derivative 
of the groundstate energy through the formula: 

1 

2 Jdttg 



|V^| 2 r e da, 



(4) 
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where in the above a is the usual surface measure induced by the canonical 
Euclidean structure W 1 , and rg is the shape derivative, i.e. the field 

rg : dn e -4 R p 

such that formally the boundary variation writes down: 

dQ.g +d e = {x + n(x)r e (x) ■ d6\x € dtlg} , 

where n(x) is the exterior normal vector at x S df^. If {ipe}g e ^ P is a smooth 
family of smooth functions such that ipg(x) = for x € d£lg, then the shape 
derivative rg can also be defined through 

rg{x)V^g{x) ■ n{x) = -Vgi>g{x) ^ x G d ^e- (5) 

Formula (5) can be proved as follows: from the Dirichlet conditions, one has 
for any small h: 

= (x + h- r e {x)n(x) + 0(|/i| 2 )) - ^(x) 
= h ■ r e (x)Vtpg(x) ■ n(x) + h ■ V e ip* e (x) + 0(h). 

Differentiability in formula (4) is a classical result of abstract analytic per- 
turbation of linear operators (see [24]), but can be proved directly with the 
variational formulation as in [19]. 

In Section 3, we introduce a standard Wiener process (Brownian motion) 

t^W u 

with some given initial distribution in Qq. The first exit time of the domain 
Qg is denoted by 

T := inf (t > 0\Wt € dflg) . (6) 

The long time probability distribution of the latter process with Feynman- 
Kac weights, and conditioned to remain in the domain Qg is denoted drjg 

r E L{W T )l T < T e- & y(Ws)ds 

1 (fidrin := lim — - — -. = r — — . (7) 



It has a probability density function given by the signed groundstate ipg: 



if Ipg dx 

^ , (8) 

/ 4>g dx 
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and the exponential rate of the evolution of the weighted extinction proba- 
bility yields the groundstate energy: 

lim -i InE ( l T<r e- & v ^) ds \ = E* g . (9) 

T— >+oo T \ ~ J 

The probabilistic interpretations (7)-(8)-(9) have some variants (see (34)- 
(35)) where a drift is added to the Wiener process and the range of the 

potential V in the Feynman-Kac weig is reduced. This leads 

to Monte-Carlo methods with some importance sampling variance reduction 
which can efficiently compute the couple [ijjt >Eq). This method has been 
widely used and studied in many fields. Special care is required to treat 
the weig ht e~fo V(W 3 )ds when 

averages are computed. For instance when 
several random processes are simulated in parallel, some re-sampling of the 
set of processes has to be carried out at regular time intervals, according 
to the weights associated with each process. We refer the reader to [21, 1, 
2] for applications in Quantum Chemistry (Diffusion Monte-Carlo or Pure 
Diffusion Monte-Carlo methods), to [16, 17] for applications in Bayesian 
statistics (Sequential Monte-Carlo methods), and to [14, 13, 15, 26] for the 
associated mathematical analysis. 

Then, we consider the weighted distribution of the Wiener process at 
the hitting time r, when the process is initially distributed according to rjg 
defined in (7). The latter distribution is denoted dfig A and reads 

f <pdf4 A := E f^(WV)l r<+00 e-/o r ^CW.)-A)«fa | Law{Wo ) = v*e) 
Jan e v ' 

E(<p(W T )l T < T<+00 e-fo(V(W s )-x)^ 



lim 



E[l T<T e- ^(V(w s )-x) ds 



(10) 



which verifies 



/ l P d Ve,\= f / tpVipQ-nda. 

Jdn e 2 (E* g — a) / r e JdQe ( ' 

Jn g 

Formula (11) holds for any A < Eg and is the grounding formula of this 
paper. Note that it can be related to the Dirichlet energy variation VgEg in 
(4) by remarking that ■ n\ = \ Vipg\. Up to our knowledge, the formula 

(11) has never been pointed out in the literature, although the sensitivity 
analysis carried out in [12] yields a similar formula but at finite time (as 
opposed to large time, which is the case here). 



4 



Approximations of V eEt with Monte-Carlo methods can then be car- 
ried out using an approximating sequence of the Dirichlet groundstate (3) 



The main limitation of computing V gEg with the Monte-Carlo technique 
suggested above is the necessity of an analytical approximation ipg of the 
groundstate ipg. Especially, the pointwise convergence of the normal deriva- 
tive VipQ.n may be hard to achieve. However, for practical situations in high 
dimension, we do not know any alternative point of view. A clear motivat- 
ing example of a high dimensional problem is the case of Fermionic systems 
where the so-called fixed node approximation is used. 

In Section 4, we introduce Fermionic groundstates (ipp,Ep) associated 
to a finite symmetry group S C 0(M rf ) of the Hamiltonian H in (1); where 
0(10^) denotes the group of isometries. S is simply the permutation goup 
of identical particles for physical systems. Fermionic groundstates are the 
solutions to the variational problem: 



ipg ^ ipg, and random samples of size approximating the probabilis- 
tic formulations (7)-(10), denoted r)f N ^°°y ^* anc j jV ~ >00 ) \- The 
following identity can then be used: 




(12) 




\ 



(13) 



Any function ip verifying the symmetry property 



VS G S, ipoS = det(S)ip 



will be called skew-symmetric, whereas any function ip verifying 



V/S 1 £ S, ipoS = i) 
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will be called symmetric. Note that existence of (ifj^Ep) follows in our 
context from the fact that H has a discrete spectrum, and is a classical re- 
sult of spectral theory for more general potential V (see references in [6]); 
but uniqueness does not hold in general. In practice, (ipp,E£.) is computed 
using a parametrization of skew-symmetric functions, and a numerical opti- 
mization procedure. This is the main problem of computational Quantum 
Chemistry, and forms a huge scientific field. We refer to [7] for a mathemat- 
ical introduction with a consequent bibliography. See also the following two 
typical papers [29, 28] involving wave function optimization using a Monte- 
Carlo method. Monte Carlo methods in computational Quantum Chemistry 
are referred to as Quantum Monte Carlo (QMC) methods. For physical sys- 
tems, the parametrization is given as a finite sum of Slater determinants 
multiplied by a strictly positive symmetric factor (called the Jastrow fac- 
tor). The result of the latter optimization relies crucially on the quality of 
the parametrization, and will be called the trial wave function, which is a 
skew-symmetric function denoted ify We will then consider a family of 
skew-symmetric functions {V'ejggjjp) with an explicit analytical expression, 
which includes ipg . The latter is used to define the nodal domains 

where: 

' M+ = jx G R d \i)\{x) > o} 

< Ng = {x € R d \ $\{x) < 0} (14) 
8Mb = {x € R d \ ift(x) = 0} . 

The fixed node approximation consists then in computing with a Monte-Carlo 
method the Dirichlet groundstate 

of the variational problem with Dirichlet conditions in J\f^ U Ma ■ We refer 
the reader to [9, 11] for historical papers on the Monte-Carlo computation 
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of fixed node groundstates. The latter variational problem reads explicitely: 



E, 



FN 



inf 



\ 



4 N H « N ) 



FN\2 



inf 



\ 



i>\dN e = 



« N ) 



(15) 



/ 



FN\2 



the same definition holding in Af g ~ by symmetry. From now on, and Mg~ 
will be called respectively the positive and negative nodal domains, dNe the 
nodal surface or nodes, and [Eg , tpg N ) the /ixed norfe groundstate elements. 
Now the key problem is the following: the energy (called Variational Monte 
Carlo, or in short VMC energy in the QMC literature) 



Ma 



4h {4) 



IN 2 



(16) 



(4) 



of the trial wave function t/jg can be minimized towards the Fermionic ground- 
state energy E* F using some optimization procedure, and the fixed node en- 
ergy Eg N < El associated with ijjg can be computed using a Monte-Carlo 
method associated with (9). However, a parameter optimizing the energy 
El of the trial wave function does not in general, for a given parametrization, 
optimize the fixed node groundstate energy E™. An open problem is now to 
develop an algorithm that can directly minimize E™. To precise this idea, 
let us consider the formulation of the exact Fermionic groundstate (i/j f , E* F ) 
as a variational problem involving the fixed node groundstate (iftg ,Ej ) 
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and the nodal surface dMe, that is to say: 

E F : = inf (^E™ , ip™ solution of (15), ipf skew-symmetric^ 

r F H(r F ) (ir) 

An approach to solve the variational problem (17) consists in the computa- 
tion of the shape derivative of the fixed node groundstate 

V e £ e FN (18) 

using a Monte-Carlo estimation based on (12). In this context, the key 
formulas (11)-(12) can be rewritten as follows, t h-> W t + denotes a Brownian 
motion in , and r + the hitting time of 8Mb . Then the probability measure 
dr]™ on Afg~ is defined by: 

E(^+)l T < r+ e-/o T ^+) rf A 

tpdif e ^ := lim — i — -. = — — (19) 

Mt T ^ Efl T<T+ e-/o T ^+)^ 



and the measure <i/i™ on dAfg for A < E™ is defined by: 

- (v(W+)-X) ds 

We will show that (12) becomes: 



E L(W+ + )l T < T+<+00 e 
/ (pdfi™ := lim -p ^r— 



(20) 



JdNg 

! r+V s ^N. n+ ^FN (21) 



(^ FN ~ A) 
2(£™ - A) 



In the above, r~g is the shape derivative of Na~ , n + is the associated exte- 
rior normal vector, and V sy V>0 N ■ n+ is the symmetrization of the normal 
groundstate gradient, that is to say: 

V S> V • n + = - (V + i) ■ n + - V"^ • n_) , (22) 
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where V +- i/' • n + (resp. X7~t/j • n_) is the exterior normal derivative in Na~ 
(resp. Nq)- By construction, r~^ is a skew-symmetric field on dJ\fg, and 
V sy ^ ■ n+ is symmetric. Our main result concerns then the link between 
(i) a symmetry breaking of the measure /u™, (ii) the fact that ipg is an 
eigenfunction, and (iii) local exrema of 9 — > E™ . The link between (i), (ii) 
and (iii) can be stated through the following equivalent assertions: 

1. The measure /x™ on dAfg is symmetric (i.e. invariant by the action of 
S). 

2. When defined on the whole space M. d , the gradient of the fixed node 
goundstate V^f N is continuous on <9AV 

3. The fixed node goundstate ipg is a skew-symmetric eigenfunction of 
H on M d . 

4. The fixed node energy variation vanishes VgE™ = for any parametriza- 
tion 9 \— > dJ\fg of the nodal surface. 

This yields a probabilistic characterization of the nodes (set of zeros) of 
skew-symmetric eigenstates of H through a symmetry argument. This is 
an original result. The practical computation of VgE™ using (21) re- 
quires a Monte-Carlo estimator of the measures (/x™, f?™) on the one hand, 
and an analytical approximation of the fixed node groundstate elements 
(?/>0 N , V sy ^™ ■ n+) on the other hand (see also (12)). A key remark is 
that the elements (ipj , V sy ip™ ■ n+) and any approximation with a skew- 
symmetric function smooth on W d , for instance with (i/jg,S/ipg ■ n + ), share 
the same symmetry properties. This suggests the following approximation 
of the energy variation VqE™ in (21): 

V^ FN ~V5™ = 2 r (E \~ll I rjv4-n + d„™ (23) 



Xv+ ^1 dr le N JdM ( 
2(E™ - A) 



Xv+ ^6 d7 l™ JdNt 
[ (H-E™) (V d 4)dn* 



[ Ve4d»Z (24) 



*3 



FN 



(25) 



In the same way as for the exact expression VgE™, the latter vanishes 

(VqE™ = 0) if 1, 2, or 3 holds. In return, if VgE™ = for any parametriza- 
tion of the nodes then 1, 2 or 3 holds. Such algorithms may be referred to as 
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Nodal Monte-Carlo. They will require variance reduction techniques exploit- 
ing the symmetry structure, in principle such that the variance of VgE™ 
scales appropriately to in the limit where //™ becomes symmetric (in other 
words, we seek for an 'asymptotically scaling variance reduction', also called 
'zero- variance estimation' in the QMC literature). Ideas are given for future 
work on this matter. 

Let us now position the content of this paper in the context of the QMC 
literature. First, we recall that efficiently computing in high dimension d 1 
the Fermionic groundstate (13) using Monte-Carlo methods is a fundamen- 
tal problem with many applications; for instance it amounts to solve the 
eigenvalue problem for excited eigenstates, where classical power methods 
fail. A general solution is known to be intractable, and is usually referred 
to as the sign problem (see Remark 4.2). This explains the necessity of the 
fixed node approximation. The issue of optimizing the nodes of the trial 
wave function ip 9 in the fixed node approximation was pointed out in [10], 
where a long discussion on the structure of Fermion nodes and appropriate 
(from this perspective) trial wave functions is provided. Yet state-of-the art 
numerical methods optimizing the trial wave function is based on either, (i) 
the minimization of the VMC energy Eg in (16), as in [29, 28]; or (ii) the 
minimization of the variance of the local energy 

E L :=V-^\)- x ^{4) (26) 

as in [29, 22]. As a consequence, an efficient method optimizing directly the 
fixed node energy E™ in (15), that is to say the nodes of the trial wave 
function i/jg, remains an unsolved problem and motivates the content of this 
paper. However, methods to approximately compute the gradient VgE™ 
were already suggested in the QMC literature in the more general context 
of the calculation of physical properties (or "forces"). The main goal is to 
compute the derivative Vi?-Ep R of the groundstate energy with respect to 
some parameter R parametrizing the Hamiltonian, typically the potential 
energy V = Vr. A classical example is the following: R is the vector of 
the nuclei- nuclei distances in molecules. The trial wave function now also 
depends on R: ipg = ify g. Well established methods are available for the 
exact and variance reduced computation of the variational energy gradient: 

Vfl4,fl (27) 

where Eg R is defined by (16). For instance, methods using coupling (cor- 
related sampling) to reduce variance were tackled in [18], and a general 
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construction of variance reduced estimators ('zero variance-zero bias' esti- 
mators) was proposed in [3]; see also some applications in [27]. The case of 
the fixed node energy gradient: 

V R E™ (28) 

where -E™ is defined by (15) is more intricate, and requires some approx- 
imation since the derivatives VriPq ^ of the fixed node groundstate remain 
unknown. Note that computing the gradients VqEq or Vgi^J^ can be seen 
as a particular case of computing respectively (27) or (28), since the nodes of 
the trial wave function ip l R e depend on R a priori . An approximate formula 
has already been proposed to compute (28), see for instance equation (54) 
in [3], (10) in [8], (10) in [5], and (14) in [4]. The terms due to the variation 
of the nodes (which amounts to evaluate V^E 1 ™) are called 'nodal Pulay 
terms', and were particularly pointed out in [5, 4]. In the above references, 
the formula used for the calculation is exactly the formula (25). However, 
the interpretation in terms of a stochastic process stopped on the nodes in 
(24)- (23), the analysis of the symmetry of the associated distribution on the 
nodes, and the suggestion of an associated Nodal Monte Carlo method are 
new results. 

The following classical textbooks are recommended: 

• about spectral theory of operators: [25]; 

• about elliptic theory of Partial Differential Equations: [20]; 

• about random processes and Feynman-Kac representations: [23]; 

• about Monte-Carlo methods in Quantum Chemistry (QMC): [21, 1]. 

2 Shape derivatives of Dirichlet groundstates 

In this section, some notations and results are recalled concerning Dirich- 
let groundstates of Schrodinger operators with a smooth potential. Then, 
formula (4) is proven formally, and references are given for rigorous proofs. 

Consider the Schrodinger operator (1) defined on M. d with a smooth po- 
tential V bounded from below. H defines a self-adjoint operator on the 
Hilbert space L 2 (M d ). For simplicity, V is assumed to go to infinity at infin- 
ity such that H has a compact resolvent, and thus a purely discrete spectrum. 
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Generalization to operators involving continuous spectrum, although of fun- 
damental importance in Quantum Chemistry, is left as technical extensions 
of the present work. Then a parametrization of smooth domains (2) is con- 
sidered for 9 £ M. p such that d£lg has a smooth boundary for any 9. One 
assumes then that there exists a set of diffeomorphisms smoothly indexed 
by 9 and such that: 

(9, x) ^ R e {x) is smooth in W x R d 
R = ld 

(29) 

Rg = Id outside some compact set 

k wew, n g = Rg(n ). 

In this setting, the shape derivative of 9 i— > Qg can be defined as the smooth 
field: 

rg : dttg -> W, 

verifying: 

Vx E dn e , r (x) = V g Rg(x) ■ n(z), (30) 
and thus locally for small h € M p : 

dQ,g + h ~ {x + h- r e (x)n(x)\x € cTi^} . 

Now classical results of spectral theory ensures that the Hamiltonian (1) 
considered in h 2 (Qg) with Dirichlet boundary condition, is self-adjoint with 
domain T>(H) C Hq(Qq), where H^Qg) is the usual Sobolev space of func- 
tion with Dirichlet conditions and square integrable first order derivatives. 
Moreover, H has a unique (up to a multiplicative constant) signed ground- 
state tjjg, solution of the variational problem (3), or equivalently solution of 
the eigenvalue problem: 

I 4>*g\dn e =° 

{ r e > o, 

where - \dCig denotes the usual trace operator on the boundary. The regularity 
of V then ensures that tjjg is smooth on £lg and that V^g • n is smooth on 
dQg. The following derivative formula can now be stated. 
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Lemma 2.1. Consider domains 9 \— > fig verifying (29). Then the Dirichlet 
energy Eg solution of the variational problem (3) is differentiate with respect 
to the parameter 9, and the variation formula (4) holds. 

Proof. The formal computation is detailed, and references are given for the 
rigorous result, which is not much more involved but less instructive from 
our point of view. Fix 9$ € W , and consider the formal chain rule: 



(VoE* g )o- 



■e 



( I m(rA 

V - 



m 2 



i 



( 



f^fl f)=8 
\ 



I 



+ 

' e=e \ 



V 6 



e=e 



Since (ipQ,Eg) is an eigenelement, it yields: 



V, 



VeE* do = 0. 



e 0J 



Then formal differentiation yields 

/ pgHW) [ v e ii>m-Eo)M) i r 9 (H-E* e )(v e r 9 ) 



V 



(Y>, 



*\2 



and since (ipg,Ea) is an eigenelement, the first term of the right hand side 
vanishes so that finally: 



(y 9 m 



[ ^* o (-lA + V-E* eo )(Ve o re ) 



(r 9o y 



Now applying Green's integration by parts two times, and remarking that 
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(H-E* o )(r eo )=0,weget 

-nVg^ da 



(VeEg) e=6o = - 



0o J 



i) 



Then (5) applied to ipQ yields the result. The rigorous proof can be made 
using a change a variable with the diffeomophism Rq, and then exploiting 
the variational formulation (see e.g. Theorem 2 in [12]). □ 



3 Probabilistic interpretations 

In this Section, the probabilistic formulations (8)-(9)-(ll) are proven and 
detailed. Associated Monte-Carlo methods are recalled with some references. 

Consider notations and assumptions of Section 2. Let t t— > Wt be a 
standard Wiener process with exit time r from defined in (6). The 
classical probabilistic interpretation of the eigenelements (ijjn ,Eg) is recalled 
in the following lemma: 

Lemma 3.1. Assume Wo is distributed according to pf^p%; where V'init € 

J Vinit \X)aX 

L 2 (Qg), and ipi n it is non vanishing in each connected component ofQg. The 
groundstate elements can be expressed through the long time behavior of the 
process with Feynman-Kac weights and conditioned by large exit times. This 
yields the formulas (7)-(8)-(9). 

Proof. The result follows from the classical representation of parabolic equa- 
tions through the Feynman-Kac formula (see [23]). Let us recall the differ- 
ent steps of the argument. First, consider ip G C£°(£Iq) a smooth solution 
(t, x) i-> u t (x) € C°°(]R + x Slg) of the parabolic problem: 

d t u t = -H(u t ) 
< u t \gu e = 
, uo = ¥ 

Then using ltd calculus, it yields for any (ft € C°°(M + x £lg): 

(f " V ~ dt ) ((/,r - 4) ( w t) e ~ J ° ViWa)dSdt + e ~ J ' ViWa)dS ^^T-t(Wt) ■ dW t , 

(31) 
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so that applying the latter computation to (f> = u between time and the 
stopping time T A r = inf (T, r) yields: 

ur(x)iphat(x)dx 

=E(u T _ TAT (W TAT )e-fo AT v(Ws)^ 

ipi nit (x)dx 

= E(ip{W T )l r > T e-fo v ^ ds 



Now denoting by (ijjg'' n , Eg' n ) n >Q the full spectrum of H normalized in L 2 (ilg) 
the above Feynman-Kac representation reads: 

E(# T )l T > T e-/. T ^)j=^e- £ »" T ^ J (32) 



™>0 / V'init 

Jn g 

Now using the fact that the groundstate has a spectral gap (Eg' 1 > El' = 
Eg), the dominant term in (32) when T —> +oo enables to verify that the 
groundstate has a sign t/jg'° > 0. Finally, taking the limit ip — > 1 by dom- 
inated convergence in the formula above, and then the leading term in the 
limit T -> +oo yields (9). (8) follows. □ 

In practice however, a diffusion solution to a stochastic differential equa- 
tion with a repulsive drift at the boundary d£l is used: 

f dX t = {^-^^{Xtjdt + dW t 

(33) 

{ ^\dn g =0, 

where ip 1 > is a smooth function strictly positive in and vanishing 
on dQg. In [6], sufficient conditions on the behavior of vp l near dflg and at 
infinity are given for (33) to be well-posed, and to verify the following variant 
of (8)-(9): 

lim _>. L = llh (34) 

lim -- InE ( e-fo E L(Xs)ds\ = E * ^b) 

T— >+oo TV / 
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where the so-called local energy is defined by (26). The proof of the latter 
probabilistic interpretation (34)-(35) is similar to the proof of Lemma 3.1, 
and is based on the mapping of the Hamiltonian H to a weighted L 2 space 
through: 

Hj = (^y l H(^-) = - - - (v'r'vv 1 -V + E L . 

Details and assumptions can be found in [6]. 

Next the probabilistic interpretation of the shape derivative given by 
formula (11) is proven. 

Proposition 3.2. Assume Wq is distributed according to JfegsiM^g where 

J Vinit \X)ax 

V'init € L 2 {VLq), and il>\mt is non vanishing in each connected component ofQg. 
Assume the boundary dVtg is smooth and uniformly Lipschitz, and consider 
the measure defined for any A < Eg by (10). Then (11) holds. 

Proof. Step 1. Let (p 6 C°°(0,g). We claim that the parabolic differential 
equation with inhomogenous Dirichlet conditions 

dtht(<p) = -(H-\) (ht(<p)) 

ht(<p)\dQ g = <p\dn e (36) 

ko((p) = if 

has a unique smooth solution for A < E* converging exponentially fast to- 
wards hoo((p) unique smooth solution of the elliptic inhomogenous Dirichlet 
problem: 

r (H-\)(h 00 (<p)) = o 

(37) 

I hoo(ip)\dn e = (p. 

This is a classical consequence of spectral theory, but we recall briefly the 
basic arguments. Existence of a smooth solution of (37) follows from the 
fact that (H — A) -1 can be extended to a bounded operator of L 2 (Oe) so 
that: 

hooitp) = (H-X)- 1 (H-X) (if) -if, 

which indeed is solution of (37). Note that in the above (H — A)~ x implicitly 
refer to the operator with homogenous Dirichlet boundary condition, so that 
(H — A)" 1 o (H — A) ^ Id when operating on test function with inhomogenous 
boundary conditions. Then the homogenous solution of (36) has been solved 
using spectral decomposition in (32), proving that such homogenous solution 
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is unique and vanishes exponentially fast when A < Eg. The latter analysis 
of the homogenous case proves uniqueness in (36)- (37), as well as exponential 
long time convergence of the time dependent equation (36). 
Step 2. We claim that for any ip £ C°°(J} d ), 

^ooM^init 

:= E L{W T )l T<+ooe - fo T (VW s )-X) ds\ (38) 



^init 

where h^tp) is the solution to the elliptic partial differential equation with 
inhomogenous Dirichlet conditions (37). Indeed, using (31) with test func- 
tion h,T-t{}P)i an d stopped at time T At = inf(T, r) yields: 

hr(<p)il>wit 

- := E (h T „ TAT &)(WrAT)e-fo AT (y(W a )-X)ds^ 



■0init 

n e 

Now the event l T=+00 has null probability if 1 1— > Wt is recurrent, and if the 
latter is transcient then using limF = +oo: 

oo 

hm E(l r=+00 e-fi AT VW)- x > da ) =0; 

so that taking the limit T —> +oo leads to (38). 
Step 3. We claim that 

hoo{fWe = — rv / ¥>V*/>e -nda. 
vie l \^e ~ A ) Jdn, 



Since ipg is the groundstate: 

re = Ef^X {H ~ X) m ' 
so that integration by parts yields: 

/ /looMV^ = T-ift ~ r / -<pVi/Jo -nda 



1 


E* e - 


A 


1 




1 e;- 


A 






E*- 


A 



-V$J • VMv?) + (V- A)Mv)$5 
-(pVipg ■ nda + 0. 
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Step 4. Now using the Markov property of Brownian motion yields: 
E(^(iy T )lT<.< + ooe-/o T (^)-^) = 

" E (E (^(^)l f<+00 e-i'o f (^(^)-^)^|^ = Wt ) l T < T e-£W W ^ ds ) , 

where ((Wt)t>o, t) is an independent copy of (Wt)t>o, t). Together with Step 
2 and Step 3 with the probabilistic interpretation of the groundstate ipQ in 
(10)-(8), it completes the proof. □ 

4 Fermion groundstates 

Consider the Schrodinger operator (1). Fermionic groundstates (tp F ,E F ) 
associated to (1) are defined with respect to a finite symmetry group S C 
0(R ) of the potential V, where 0(M d ) denotes the group of isometries of 
M d . A symmetry group of V is defined by the property 

VS* € S V o s = V. 

Fermion systems appears in the special case when S contains symmetries 
with odd parity: 

det(S) = -1, 

so that H can be defined (since the Laplacian operator commutes with isome- 
tries) to operate on the Hilbert space of skew-symmetric function: 

U = {V> € h 2 (R d ) | VS G S, ^oS = det(S)^} . 

"Fermion" groundstates are then the solutions to the variational problem 
(13). 

Example 4.1 (Fermions). In the case of iV physical quantum particles of 
Fermionic type in dimension 3, one has R d = R 3N , and a potential of the 
form: 

N 

V(x) := ^2 V ext (xi) + ^2 Vint{%i ~ Xj), 

i=l l<i<j<N 

where x = (x±, . . . ,X]y) are the 3 dimensional coordinates of the particles, 
Vext is an exterior smooth potential that goes to infinity at infinity, and V[ nt 
a smooth interaction potential that vanishes at infinity. Then the discrete 
symmetry group is the permutation group Sjy of physical particles. See [6] 
for the mathematical analyis of Quantum Monte-Carlo (QMC) methods in 
this context. 
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Since in this paper, we restrict to smooth operators with compact re- 
solvent, the existence of (ip* F ,E* F ) follows directly from the fact that the 
spectrum is discrete, and smoothness from the fact ip* F satisfies an eigen- 
value problem of a smooth elliptic operator (see [20]). 

Remark 4.2 (The sign problem). Computing directly (ip F , E F ) using Monte- 
Carlo methods is an untractable problem known as the sign problem. The 
latter can be summarized as follows. Since H leaves invariant 7i, eigenfunc- 
tions of H in M rf are either symmetric or skew-symmetric functions. Thus 
computing (32) with: 

• a skew-symmetric test function ip l , 

• a non-symmetric positive initial condition t/'init > 0, 
yields: 

/ ^"^it , 

ip*,«en / Vinit Rd 



and in principle if <p is symmetric: 

e ( K y{w T )i> l {W T )<r ! ° v{Ws)ds ) J Kd ^f^V 



Unfortunately, this computation relies crucially on the rate of vanishing of 
the normalisation which is due to the non-symmetry of the initial condition 
^init > 0. Since stochastic processes used to compute such averages quickly 
forget their initial condition and have a symmetric dynamics in the full state 
space one is compelled to compute ratios of vanishing averages (of the 
type |j), with Monte-Carlo estimators having a non- vanishing statistical vari- 
ance. This leads to infinite variance when T is large. This forms the sign 
problem. This problem appears more generally when trying to solve higher 
eigenvalue problems with Monte-Carlo methods. Although, there is proba- 
bly no general solution, solving the sign problem for particular situations in 
high dimension would be considered as a major breakthrough. 

In practice, (ip F ,E F ) can be approximated using a hybrid methodology 
in two steps. First, a trial wave function is obtained using an analytical 
parametrization of H, usually of the form: 

<e ■= « W (39) 
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In the above, J a > is a strictly positive symmetric part called the Jas- 
trow factor and parametrized by a G M n ; and ^ kew is the skew-symmetric 
part parametrized by 9 € MP. The most convenient numerical optimization 
method is then used to solve (13) in the space of parameters formed by 
(a, 9), and the solution of the optimization procedure is denoted with the 
parameters (ao,#o)- This yields the set of functions {V'ejggRp defining the 
nodal domains through 

4 := J Q0 ^ kcw 

Example 4.3 (Fermions). In the case of N physical particles of Fermionic 
type, Vf cw is built using a sum of Slater determinant, that is to say a sum 
of functions of the form: 



det 



bj\ x Wi,j=i...N 



where ((f) j) ■ 1 N are N smooth functions of Mr. 

The link with Section 2 and 3 is made by posing: 

fie := Mq UjV r 

The Fixed Node Approximation consists in computing with a Monte-Carlo 
method the solution (ipj , E™) = (i!jq,Eq) of the variational problem with 
Dirichlet conditions (15). Such a computation is made using the probabilistic 
interpretations (8)-(9), or more usually in practice using the variant (34)- 
(35). This method is known under the DMC acronym (Diffusion Monte- 
Carlo) in Computational Chemistry, and has been widely studied, see for 
instance [21, 2]. 

As explained in the introduction, the key problem is that the nodal sur- 
face dJ\fg = (tpg)~ l (0) may be different from (^) _1 (0). The open question 
is thus now to carrry out numerical methods associated to the variational 
problem (17). In this context, a direct minimization of E™ in (17) requires 
the computation of: 

VqEq 

using formula (4) and the probabilistic interpretation (11). 



5 Fixed node and symmetry breaking 

In this section, we consider the context of Section 4, and we assume that 
the nodal surface defined in (14) is a smooth manifold, a sufficient condition 
being: 

Vx G dMe, V^J(x) ± 0. (40) 
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Moereover, we assume that the assumptions of Section 2 on the mapping 9 
Mg apply, namely that there is a diffeomorphic mapping 9 i-> Rg associated 
to the domain Me such that (29) holds. 

A fundamental remark concerns the symmetry of the normal derivative 
of the fixed node groundstate S/ + tp™ .n + , (where V + and n + refer to the 
domain Mq~ , and will be defined below) , or equivalently the measure fig A 
in (11) defined on dMg for stochastic processes evolving in M^ . The latter 
indeed presents a symmetry breaking, in the sense that they are only invariant 
by the action of the symmetry sub-group 

S + := {S 6 S | det(S') = 1} 

on dMg. Before going further, we will precise notations in appropriate defi- 
nitions and lemmas. 

Lemma 5.1. S is a symmetry group of the nodal surface dMg in the sense 
that any space transformation S € S verify: 

S(8Mg) = dMg. 

Proof. By skew-symmetry, it yields for any x £ M. d 

4(s(x)) = -4( X ) 

so that ipg(x) = is equivalent to i^ l d {S{x)) =0. □ 

So consider now ip a skew-symmetric function in M rf with Dirichlet bound- 
ary conditions on the nodes ip\dAf e = 0, and such that the restrictions 
ip\j\f+ '■ M~g^ — > R, on the one hand and ip\j^- '■ Mg — > R, on the other 
hand are smooth on the associated closed domains. Remark that using 
skew-symmetry, for any S € S verifying detS" = — 1, tp\j^- is the image of 
ip\j\f+ through the transformation: 

Two traces of on dMg can then be defined depending if the reference 
domain is M^ or Mn~ . The definition of associated symmetric and skew- 
symmetric traces then follows. 

Definition 5.2. Let ip a skew-symmetric function in R rf with smooth restric- 
tions and Dirichlet boundary conditions on Mg and Mg . n + denotes the 
exterior normal of dM g + , n_ denotes the exterior normal of dM d ~ , so that: 

n + = — n_. (41) 
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V + ip ■ n + denotes the exterior normal derivative of tp in Ma~ , and V _- i/> • n_ 
the exterior normal derivative of tp in M e ~ . Then the skew-symmetrization 
of the normal derivative is defined by 

V sk V> • n = - (V + V> • n+ + V - ^ • ra_) , (42) 

and the plain symmetrization with respect to n + is defined by (22) . 

Note that tp 6 C 1 (R d ) if and only if V sk ^ • n = on dM e , so that V sk ^ • n 
can be seen as the gradient discontinuity of tp at dMe- One can now precise 
the idea of symmetry breaking on dMe- 

Lemma 5.3. Let tp a skew- symmetric function in R rf with smooth restric- 
tions on Me and Me , and Dirichlet boundary conditions on dMe- Then 
V sk '^/ , • n is skew-symmetric, and V sy tp ■ n + is symmetric, under the action 
of S on dMe- If <P G C£°(R d ), one has the integration by parts formula: 

If V(p-V^= f <pV sk tP-nda- [ (^(VO- (43) 
z JR d JdN g JR d z 

Moreover, V tp ■ n = on dMe if an d only if • n + or \/~tp ■ n_ are 

symmetric. In the opposite case, V + tp ■ n + and V~tp " n - are invariant only 
under the action of the special sub-group S + and a "symmetry breaking" 
occurs. 

Proof. By skew-symmetry of tp, one has for any 5 € S in Mg~ UM g ~: 

SVipoS = det(5) Vtp, 

and then on dMe, 

S n + o 5 = det(S) n+ 
5n_o5 = det(S) n_, 
so that since S T S = Id, it yields on dMe for any S € S : 

• ?i+ if det(S) = 1 
V~tp • n- if det(5) = 1 (44) 
-V-^-n_ if det(S) = -1. 

This yields the symmetry properties of Vtp sk -n and Vtp sy -n+ with respect 
to 5, and of V + tp ■ n + and V~tp • n_ with respect to S + . 



< V^V • n_ o S = 
, V + V> • n + o S = 
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The integration by parts formula is obtained by applying seperately on 
Af g + and M d ~ the classical Green's identity. 

Finally from the symmetry properties (44), V~ip ■ n_ = — V + tp ■ n + on 
dNg if and only if V + ip ■ n+ o S = • n + or V - ^ -Ti—oS= V - ^ • n_ for 
any S G S with det(S) = -1. □ 

One can now apply these remarks to the fixed node groundstate ip™ 
solution of (15). 

Lemma 5.4. Let ip™ be the solution of (15) with a smooth boundary dJ\fg. 
Then the following assertions are equivalent: 



1. ip™ is an eigenfunction of H in L 2 

■ n, as defined by (42), vanishes on dNg. 

3. V + ^™ • n + or X7~t/j™ ■ ra_ are symmetric under the action of S on 
dN e . 

Proof. The third point and the second point are equivalent by Lemma 5.3. 
Now, if ip € C^°(R d ), integration by parts (43) yields: 



so that ip™ € L 2 (IR d ) is an eigenfunction if and only if V sk, i/'™ • n = on 
dAfg. □ 

In the present context, the formula of the shape derivative of the Dirichlct 
groundstate can be written with symmetrized normal derivatives. We denote 
the shape derivative rt (resp. r^) of the nodes as defined in (30), with 
respect to the normal n + (resp. ra_), so that 

r+ = -r~. 

Proposition 5.5. The shape derivatives of the fixed node groundstate ifj™ 
with respect to the nodal parameter 9 reads 

V e £™ = — * / (VV™ • n) (V s ^ e ™ • n+) r+ da. (45) 

/ m 2jm ° 

Jn+ 

Moreover, VgE™ = for any parametrization 9 i— > dflg of the nodes if and 
only if the fixed node groundstate -0g N is an eigenstate of H in M. d . 



23 



Proof. The shape derivative formula (4) can be decomposed as the sum of 
the part due to Af^~ and to M e ~: 

VoEj" = * ( [ (|V^™| 2 r+ + IV-^Tfrj) da) . 

2 / (r e ) 2 yJme J 

Since ipg = on cWg, then |V + V FN | = |^ 7+ V' FN an d symmetrically 

|V-^ N | = |V"^ FN ■ nJ. We get: 



4r+ (V sk ^™ • n) (V s >f N • n + ) = r+ (|V+^ FN | 2 - |V^ FN 



2^ 



= ^|V + ^ FN | 2 + r-|V-^ FN | 2 , 

where in the last line we have used rt = —fg. Then (45) follows. As a 
consequence, \7qE™ = if and only if V + ^ FN • n + = V~^ FN • n_, which is 
equivalent by Lemma 5.4 to the fact that ip™ is an eigenfunction. □ 

We can now state the main result of this paper, which consists in a char- 
acterization of the nodes of eigenstates through the symmetry of a random 
stopped process, and suggests a method to evaluate the shape derivative of 
the fixed node energy V$Eg . 

Theorem 5.6. Consider a Brownian motion t H > in J\fg~, and r + the 
hitting time of 8Nq. Consider the measure /x™ on 8Nq defined for any 
A < E™ by (20) . Then the fixed node groundstate ^ FN is an eigenfunction 
of H in M. d , if and only if ^™ is invariant under the symmetry group S. 
Moreover, the shape derivative of the fixed node energy E™ is given by (21). 

Proof. The proof consists in a direct application of Lemma 3.1 and Propo- 
sition 3.2. Indeed, (11) yields the identity of measures on dj\fg: 



FN V + V> FN • n + da 



2{E™-\) « N ) 2 
V+ 



and Lemma 5.4 enables to characterize eigenfunctions of H from the sym- 
metry of f?J N,A . Next, since on 8J\fg is skew-symmetric, and the integral 
on 8Mb of integrable skew-symmetric functions vanishes, it yields: 

f r+V sy V> FN • n+V+V> FN • n+ da = I r+V sy V> FN ■ n + V sk ^ N ■ n da. 

Then (45) with (34) yield the result (21). □ 
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6 Comments on Monte- Carlo methods 



Theorem 5.6 suggests a Monte-Carlo general strategy to approximate the 
fixed node energy variation (21) by using the probabilistic interpretation 
(20)- (19). As already commented in the introduction, computing (20) with 
Monte-Carlo methods is a well known topic, and computing (19) can be 
carried out using independent stopped processes. In fact, the most straight- 
forward limitation of the method consists in approximating on the nodes 
dMg the symmetrized gradient 

V s ^ e FN ~ We 

by a sequence ipg n ~ 5,+ °°> -0™ of smooth approximating functions in M. d , usu- 
ally converging in the sense of the energy norm. Indeed, there is no reason a 
priori that the gradient Vipg converges in a pointwise sense towards V sy 0™, 
which is necessary when approximating the measure r/™ ,X using a Monte- 
Carlo sample. However, in the context of optimization, it is fundamental to 
remark that the stationary nodal surfaces solution of VqE™ = for any 
shape derivative are not modified when approximating V sy 0™ ■ n + by any 
symmetric field V^g ■ n + . This is the meaning of the following proposition. 

Proposition 6.1. Consider any smooth skew- symmetric function tpg in M. d 

with zeros given by dM$. Then the approximation VgE™ given by (23) 
yields (24) -(25). Moreover the fixed point equation 

V^™ = 

holds for any shape variation r^ of the nodes dJ\fg if and only if ip™ is an 
eigenfunction of the Hamiltonian H inM. d . 

Proof. Since dJ\fg is defined implicitely by ipg = 0, deriving ipQ +h (x+r~^ +h n + ) = 
with respect to h yields on dAfg: 

Ve4 + r+VV4 • n + = V e 4 - r+ | V^| = 0, 

where in the last line we have used the fact that • n+ = — | | on 8Mb 
since M e + = {iGK d ij}\{x) > 0}. This gives (24). Next (11) yields: 



9AT+ r)( ZTiFN 



2{E^-\) 



— / V 4^e N ■ n+ da, 

v) / ^FN J 8M e 
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so that integration by parts in gives (25). 

Finally remark that • n_ is symmetric and that the shape derivative 
r~a is spanning all skew-symmetric field, so that the fixed point equation is 
equivalent to the fact that fi™ is symmetric. Then Lemma 5.4 enables to 
conclude. □ 

Eventually the following algorithm (without details) is suggested to es- 
timate V e E^ N . It might be referred to as Nodal Monte-Calo (NMC). 

Algorithm 6.2 (NMC). Consider a parametrization of wave functions of 
the form (39). Then the following steps are suggested: 

1. Optimize the symmetric "Jastrow" factor J a according to the varia- 
tional problem with fixed node (15), and denote ifi the obtained trial 
function. 

2. Generate a sample according to the probability distribution r/™ in Afg~ 
defined in (20). Use for instance a long time trajectory of a drifted 
stochastic process of the type (33) with Feynman-Kac weights as in 
formula (34). 

3. Use the latter sample as initial conditions; and sample independent 
Brownian motions. Stop them when they hit the nodes dj\fg. Then 
compute the measure /U™ according to (20). 

4- Estimate the variations of the fixed node energy V e E™ with (24). 

As a conclusion, we mention two possible strategies of variance reduction 
that will be necessary for efficient computations. Their development is left 
for future work. 

• First, a drifted diffusion instead of plain Brownian motion may be 
used in (20) to reduce the variance caused by the exponential weights. 
However, the classical stochastic differential equation (33) where if) 1 
vanishes at the nodes cannot be used. Indeed, the repulsive drift at 
the nodes prevent the process from hitting the latter, and (20) no 
longer holds. Instead of a skew-symmetric guiding function we 
propose to use as a drift a symmetric function i/jr, for instance an 
approximation of the Bosonic groundstate tp^ > of H in R rf solution 
of the eigenvalue problem 
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This may remove the variance of the exponential Feynman-Kac weights, 
while letting the walkers hit the nodal surface. Such variance reduc- 
tions are referred to as importance sampling methods. Note that in 
known algorithms (e.g. DMC), such guided walkers will very quickly 
hit the nodal surface, causing high variance branching. However, in 
the proposed method, the algorithm is stopped when walkers have hit 
the nodal surface once, and no branching is performed, avoiding such 
kind of variance instability. 

• Second, a key point would be to develop a coupling (in the probabilis- 
tic sense) between the random process used to compute fig in (20) 
and another random process stopped on the nodes dftfg, with a distri- 
bution denoted JlJ . The goal is to design /I™ so that the following 
two features hold: first computing the average (24) with Jl™ always 
yields 0; second the coupling is perfect in the special case when fi™ 
is symmetric (that is to say the Monte-Carlo method computing fi™ 
and Ji™ has to be the same with the same random numbers). In the 
physics terminology, this would yield a zero variance estimator, which 
we prefer to call asymptotic variance reduction, in the sense that thanks 
to the variance reduction, the variance of the estimator scales with the 
quantity to be computed when the latter becomes small. 
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